Free energies of ferroelectric crystals from a microscopic approach 
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The free energy of barium titanate is computed around the Curie temperature as a function 
of polarization P from the first-principles derived effective hamiltonian of Zhong, Vanderbilt and 
Rabe [Phys. Rev. Lett. 73, 1861 (1994)], through Molecular Dynamics simulations coupled to 
the method of the Thermodynamic Integration. The algorithms used to fix the temperature (Nose- 
Hoover) and/or the pressure/stress (Parrinello-Rahman), combined with fixed-polarization molec- 
ular dynamics, allow to compute a Helmholtz free energy (fixed volume/strain) or a Gibbs free 
energy (fixed pressure/stress). The main feature of this approach is to calculate the gradient of 
the free energy in the 3-D space {Px,Py,Pz) from the thermal averages of the forces acting on the 
local modes, that are obtained by Molecular Dynamics under the constraint of fixed P. This work 
extends the method presented in [Phys. Rev. B 79 (2009), 064101] to the calculation of the Gibbs 
free energy and presents new features about the computation of the free energy of ferroelectric 
crystals from a microscopic approach. A careful analysis of the states of constrained polarization 
is performed at T=280 K (« 15-17 K below Tc) especially at low order parameter. These states 
are found reasonably homogeneous for small supercell size (L=12 and L=16), until inhomogeneous 
states are observed at low order parameter for large supercells (L=20). The effect of this evolution 
towards multidomain configurations on the mean force and free energy curves is shown. However, for 
reasonable supercell sizes (L=12), the free energy curves obtained are in very good agreement with 
phenomenological Landau potentials of the litterature and the states of constrained polarization 
are homogeneous. Moreover, the free energy obtained is quite insensitive to the supercell size from 
L=12 to L=16 at T=280 K, suggesting that interfacial contributions, if any, are negligible at these 
sizes around Tc. The method allows a numerical estimation of the free energy barrier separating the 
paraelectric from the ferroelectric phase at Tc (AG ^ 0.012-0.015 meV/5-atom cell). However, our 
tests evidence phase separation at low temperature and low order parameter, in agreement with the 
results of Troster et al [Phys. Rev. B 72 (2005), 094103]. Finally, the natural decomposition of the 
forces into onsite, short-range, dipole-dipole and elastic-local mode interaction allows to make the 
same decomposition of the free energy. Some parts of this decomposition can be directly calculated 
from the coefficients of the Effective Hamiltonian. 
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I. INTRODUCTION 



Two main theoretical approaches are currently used in the description of ferroelectricity in perovskite oxides. The 
first one, that we will call the microscopic approach, is based on first-principles calculations and effective hamiltoni- 
ans. These eflFective hamiltonians are simplifications of the energy landscape of ferroelectric solids in terms of relevant 
degrees of freedom (local modes, displacement modes, homogeneous strain and, for more complex perovskites, anti- 
ferrodistortive modes!) . They have even been recently extended to describe multiferroics by accounting for magnetic 
moments as degrees of freedom^.. Their mathematical form contains adjustable parameters that are usually derived 
from density-functional calculations. The effective hamiltonians, once constructed, can be solved either within Monte 
Carlo (MC) or Molecular Dynamics (MD) simulations. The physical quantities (polarization, strain...) are obtained 
as thermal averages over the equilibrium trajectories constructed from these two methods, that should provide equiv- 
alent sampHng of phase space (MD) or configuration space (MC). These effective hamiltonians contain most of the 
thermodynamics of ferroelectric crystals and, in most cases, describe very well their evolution with temperature and 
pressure/stress with or without external electric field. They have been extended to treat low-dimensional systems 
(thin films, nanowires^, dots^i^) under various mechanical and electrical boundary conditions. The first ab initio- 
derived effective hamiltonian has been constructed by Zhong, Vanderbilt and Rab^^ and successfully appHed to 
barium titanate. 

The second approach, that we will call the macroscopic approach, or phenomenological approach, is based on the 
Landau theory of phase transitions. This theory uses as central concept an incomplete thermodynamic potential 
(that can be a Helmholtz potential F or a Gibbs potential G, or the thermodynamic potential of any other ensemble). 
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commonly called Landau free energy, which is a function of the polarization P (order parameter) and can be defined 
from the restriction of the true thermodynamic potential to the microscopic states having a given polarization P— . 
In the canonic ensemble, it writes: 

F{N,V,T;P) = -kBTlnZ{N,V,T;P), 

with 

Z[N,V,T; P) = ^e-''/'=«^, 
i/p 

in which the incomplete canonical partition function defines a Helmholtz free energy F . The sum in this expression 
is over all the microstates i having a polarization P. Within this definition, the incomplete free energy F{P) appears 
as directly related to the density of probability of the order parameter n(-P) byi^: 

n(P) = n(o)e-[^(^)~^(°)i/'=«^ (1) 

The same definitions in the isothermal-isobaric ensemble allow to define an incomplete Gibbs free energy: 

G{N, P, T; P) = -ksTlnAiN, P, T; P), 

with 

A(iV, P, T;P)^J2 e-(^'+^^-)/'=-^, 
i/P 

the (incomplete) isothermal-isobaric partition function. 

The Landau approach describes quite well the thermodynamics of many ferroelectrics. It was applied in a pioneering 
work by Devonshire on bulk barium titanate a long time ago^ i^°i^^ and has been applied since then on many other 
ferroelectric materials and in many other forms than bulki^. In Landau theory, the thermodynamic potential is 
expanded in power series of the polarization P (the order parameter in ferroelectric crystals), usually up to 6*'* or 8*'' 
order. The quadratic coefficient is assumed to vary linearly with temperature while the higher-order coefficients are 
usually assumed as constants (more complex potentials with temperature-dependent high-order coefficients have been 
recently propose d^^i^^ ) . These free energies are used to model the thermodynamics of the system considered around 
the critical temperature and even further, even though the Landau approach is supposed to break down near the 
critical temperature. Since such polynomial developments are by nature approximate, it appears legitim to try to use 
the microscopic models introduced hereabove to get insight into the Landau potentials by using standard numerical 
simulation tools. 

However, the connection between the two approaches is not simple. Indeed, on the one hand, the MC or MD 
methods are unable to give directly access to the thermodynamic potential, that can not be computed as the thermal 
average of some microscopic quantity. Indirect methods can nevertheless be employed: using the fact that the Landau 
free energy is directly related to the density of probability of the order parameter, it may be possible to compute 
the free energy from the probability distribution obtained from a long enough MC or MD runM^. However, if such 
calculations are possible and efficient in the simple case of the (j)'^ modeU^, it is not the case of ferroelectric crystals 
described by the effective hamiltonians introduced above, for which the high-free energy parts (for instance the part 
of F at low order parameter in the ferroelectric phase) are not sampled correctly (or even not at all) , making the free 
energy accessible only in the vicinity of the equihbrium value of the order parameteriA. 

A clever alternative method was proposed by Iniguez et to calculate from MC simulations the coefficients of the 
Landau free energy expansion, that consisted in applying an external electric field in order to displace the minimum of 
the free energy in the {Px, Py,Pz) space. Recently, we proposed to use the method of the thermodynamic integration^^, 
coupled to MD simulations, to access numerically to the free energy as a function of polarization (calculated as a 
Helmholtz free energy, under constant strain tensor) . This method allows a correct sampling of the high-free energy 
regions (where the "reaction coordinate" has a low probability) since it consists of MD runs under the constraint of 
fixed reaction coordinate. Thermodynamic integration is a technique commonly applied in computational chemistry 
to obtain free energy profiles along a reaction coordinate. In the present case, the order parameter (polarization) is 
taken as the reaction coordinate. It has also already been appHed in the framework of Landau theory to study the 
cubic-tetragonal phase transition in Zr02 by Fabris et a^. Other techniques such as the Wang-Landau method can 
also be usedi^. 
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On the other hand, the computation of Landau potentials is dehcate for various theoretical reasons. Indeed, fixing 
the order parameter to a given value {i.e. sampling the phase space under the constraint of fixed order parameter) 
can raise unexpected problems, especially at low value, where the system might evolve to inhomogeneous states and 
separate into "domains" below Tc. In the case where this phase separation occurs, changing the order parameter from 
zero to higher values simply results, at least at low order parameter, in some domain wall motion (or more complex 
interface motions and reorganizations) , and the computed free energy does probably not correspond to the free energy 
expected as the central concept of Landau theory, in which sufficiently homogeneous states are assumed all along the 
free energy curve. Troster et aM- have studied this phenomenon in details in the framewotk of the so-called 0** model 
through Wang-Landau computations of the free energy. Moreover, if surface or interfacial contributions enter the 
computation of the thermodynamic potential, the calculated object might not be a volume quantity and not have the 
extensivity required by such a thermodynamic function. 

Essentially, such problems are related to the basic fact that defining properly the polarization as a continuous 
physical quantity requires the choice of a spatial averaging lengtblS (as it is the case, in fact, for all macroscopic 
quantities used in electrodynamics) . This averaging length scale L should be chosen below the correlation length ^ so 
that the local order parameters (electric dipoles in our case) remain correlated within the averaging volume and that 
each dipole only fiuctuates around the averaged value all over the averaging volume. Conversely, it should be also 
chosen much larger than the lattice parameter a (elementary distance between first neighbor local order parameters) 
so that the order parameter as a slowly varying macroscopic quantity keeps a sense. For these reasons, it is admitted 
that the Landau free energy should be defined with respect to a given averaging length L (coarse-graining length)^^, 
that should fullfil the condition a << L << ^. 

The correlation length is expected to decrease with temperature above Tc and to increase with it below Tc. However, 
the notion of correlation volume/length is quite complex in ferroelectric systemsi^ii^, due to the peculiar nature of 
the dipole-dipole interaction, that is long-ranged and leads to strongly anisotropic correlations. As a consequence, it 
is the shape rather than the volume that defines a region within which dipolar correlations can produce ferroelectric 
order—. Lines^^ pointed out more than thirty years ago the difficulty to define properly a correlation length in a 
ferroelectric. Anyway, we will assume in this paper that the notion of correlation length/volume keeps a sense in 
ferroelectric solids, that will help us to understand from a qualitative point of view the free energies computed with 
various supercell sizes and for various temperatures. 

From a numerical point of view, we use as averaging volume the supercell used for simulating the bulk system within 
periodic boundary conditions (L^). Computing a relevant free energy as a function of order parameter would therefore 
require the supercell size to be lower than ^. If such a condition is not fullfilled, inhomogeneous configurations are 
expexted to appear at low order parameter and low temperature. Since ^ decreases rapidly when the temperature 
decreases below Tc, and that, of course, a too small supercell can not be used to prevent from too large finite-size 
effects, it suggests that computing a relevant free energy within the present method is possible only above a certain 
temperature in the ferroelectric phase. However, the identification of such a free energy with the Landau free energy 
is still a controversial issue^^. In the same manner, at a given temperature, phase separation is expected to occur at 
low order parameter for large enough supercells of size L > ^. We add that even when L < ^ (around the critical 
temperature for instance, where ^ can reach high values), the computation suffers from finite-size effects, unavoidable 
in that case within periodic boundary conditions. 

In the present work, we present some features that may help to get insight into Landau potentials of ferroelectric 
materials from a microscopic approach, at least around the Curie temperature. We extend the approach of Ref. [13 
to MD simulations under constant pressure, that allow to compute a Gibbs free energy as a function of polarization. 
Once applied to BaTiOa, the free energy curves obtained are very similar to what is expected from Landau theory of 
first-order phase transitions and, after fitting these curves by 8"* order polynomial functions, an excellent agreement 
is found with availables phenomenological Landau potentials of the litterature. 

We also conduct a very careful analysis of the states of constrained order parameter for the lowest temperature 
studied (« Tc - 15 K). They are found to be more and more inhomogeneous as the supercell size increases. For 
the largest supercell (L=20), phase separation starts to occur at this temperature. Below this size (between L=12 
and L=20), the mean force and free energy curves obtained are independent on the supercell size, showing that the 
computed free energy is free of interfacial contribution and is actually a volume quantity. For a reasonable supercell 
size (L=12), the macroscopic states of constrained polarization are shown to be homogeneous (this is carefully tested 
by increasing the number of time steps). However, we show that below a given temperature (« 240-250 K) the use of 
a 12 X 12 X 12 supercell leads to inhomogeneous states at low order parameter, as expected. 

We also discuss the contributions to the free energy of a ferroelectric crystal in terms of the four contributions 
(onsite, short-range, dipole-dipole and elastic-local mode interaction) that constitute the effective hamiltonian. We 
show that two of them are independent on the temperature when expressed as functions of the mean local mode u 
(dipole) instead of polarization P. 
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II. THEORETICAL AND COMPUTATIONAL DETAILS 



A. Hamiltonian 



We have used the Effective Hamiltonian of Zhong, Vanderbilt and Rabe^*^, devoted to BaTiOa. This hamiltonian 
uses as degrees of freedom the "local modes" Ui (local order parameters), that roughly represent the local dipoles that 
do exist instantaneously in each unit cell, the (mechanical) displacement modes Vi, from which the inhomogeneous 
strain tensor {??/(«)} is constructed for each unit cell i, and the homogeneous strain tensor {vi^}- The strain tensor 
at cell i is thus given by 'rii{i) = rj^ + 77/ (i). 

This hamiltonian H'^^^ {{ui} , {rii{i)}) consists of several terms, among which a long-range interaction is included to 
model the dipole-dipole interaction (this term is absolutely necessary to provide a reaHstic simulation of a ferroelectric 
system): 

H^ff{m , {mm) = E^-'^im) + e''p\{u,}) + E^^^^^iu,}) + i?^'-({,7K0}) + i?^"*({«j , {mm) (2) 

The first term E'^'^^^ — E{ui) is local and consists of a fourth-order polynomial function in the uf. 



E{ui) = K2 \\ui\\ + a \\u^\\ + -l{u\^u\ + u\u\^ + u-^u-J (3) 

The second term E^^^ is crucial to the description of ferroelectric systems: it is the (long-range) dipole-dipole 
interaction between the local modes: 

^dvl ^^Y^ - S{Rij.Ui){Rij.Uj) 

with Rij = Rij/Rij, Rij being the lattice vector joining cell i to cell j. Z* is the effective charge associated to the 
local modes and Coo the electronic dielectric constant. This long-range interaction can be expressed as: 

^dpl ^ Qij^af3U^aUj[3, (5) 

the Q matrix being calculated by using the well-known Ewald summation method. 

The third term describes short-range interactions between neighboring local modes (up to the third neighbor): 

i?^'""^*-^EE%-/5"-"^/3 (6) 

The fourth term is the elastic energy associated to the strain. This term is calculated from the three elastic constants 
Bii, Bi2 and B44 of the parent cubic phase. 

Finally the fifth term describes the coupling between the local modes and the strain. In ferroelectric systems, this 
coupling is at the origin of the so-called electrostrictive effects and is of primary importance. It is supposed to be 
local and takes the following form: 

= ^ E E Blc.0mm^aUJl3 (7) 
i la,0 

All the parameters of this hamiltonian are determined from first-principles calculations of well-chosen 
configurations^ in the framework of the Local Density Approximation (LDA) to Density Functional Theory (DFT). 
The reader is refered to Ref.Qfor their precise value. This effective hamiltonian describes very well the thermodynam- 
ics of barium titanate, especially its complex sequence of phase transitions (rhombohedral - orthorhombic - tetragonal 
- cubic), with a Curie temperature however a little too low (around 300 K). 



B. Numerical method 



We perform MD simulations based on this hamiltonia n^^i^° . According to the algorithm used, we can perform 
either fixed temperature simulations (Nose-Hoover algorith m^^i^^ ) , fixed stress tensor simulations (Parrinello- Rahman 
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algorithmSi) or both. The combination of the two allows to reproduce with a very good accuracy the sequence of 
phase transitions and temperature evolution of strain and polarization of barium titanateiS,, as found from MC 
simulations^. The equations of motion in the Parrinello-Rahman scheme are recalled hereafter: 

respectively for the local modes {ui — H{t).u'j) and diplacement modes {vi — H{t).v'j). mim and m^sp are respectively 
the masses associated to the local modes and to the displacements modes. H(t) (the 3x3 matrix formed by the 
components of the three vectors a, b and c that define the supercell) evolves according to: 

^^rd'^H , . . . 

in which is a "mass" associated to the dynamics of the supercell vectors a, b and c, which has to be correctly 
chosen. = dVl/dHij (fi is the volume of the supercell) and G{t) =* H[t).H{t). a and oq are respectively the 
instantaneous stress tensor and the desired stress tensor. 

In the preliminary work of Ref . l20i. we have performed MD simulations under fixed polarizatiorM- . These simulations 
are achieved by adding external forces in the time-evolution equations of motion of the local modes (see Ref. [^^ for a 
precise explanation): 

with C,a = ^I^Tlii fi^- '^i,a IS the a-component of the z*'* local mode (which has the dimension of a displacement), 
mim is the mass associated to the local modes, the a-component of the force acting on the i*'' local mode and 

N the number of 5-atom unit cells in the supercell used for the simulation. The quantity C = (CxiCy^Cz), calculated 
in Ref. [IS as a Lagrange multiplier, is formally equivalent to an external electric field that would be applied to the 
system and would maintain invariant the polarization. The total instantaneous external force applied to the system 
to maintain the polarization fixed is thus — ff^ . This Molecular Dynamics under constraint is assumed to provide 
a correct sampling of the subspace of phase space corresponding to a fixed value of the order parameter, defined by 
an average over the whole simulation box. 

Based on such simulations, we have applied the method of the thermodynamic integration to compute the difference 
of free energy AF{N, {77} ,T;u) = F{N, {77} , T; u)—F{N, {ry} ,T;u — 0) between u — and u. We use for convenience 
in the following the average local mode It = l/Nj^i'^i as order parameter instead of the polarization P (P = NZ*/flu, 
Z* being the effective charge associated with the local modes and ft is the volume of the supercell) . 

This difference is obtained as the integral over a continuous path in the 3-D space {ux,Uy,Uz) of minus the thermal 
average of the total forces acting on the local modes, this thermal average being computed under fixed volume (more 
precisely under fixed strain tensor), fixed temperature and fixed mean local mode: 

F{N, {r^} ,T-u)- Fo{N, {r;} , T) = - / J] (/|™) .du' , (8) 

i 

The subscript in this formula and in the following is used to denote the quantities at u = 0, and the integral is 
performed over a continuous path joi ning m = to w in the 3-D space {ux,Uy,Uz). A complete and rigorous proof 
of this equation can be found in Ref. |20| in terms of partition functions, but basically it can be recovered from the 
principles of thermodynamics: indeed under fixed volume V and fixed temperature T, the infinitesimal variation of 
the free energy dF when the polar displacements evolve from m to u -I- c?m is equal to 6W' , the work of the external 

forces (other than pressure forces), i.e. in our case the work of the external force ~N / A = — J2i (/i" 



dF{u)^SW' = -J2{fl"^} du (9) 

"^-^ \ / N,{ri},T;u 
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which writes in a local form: 



E(/]"%,^,,^, = -v.^(^^^M^^;") (10) 

i 

and provides by integration Eq. [8l The free energy is thus computed as a potential of mean force (PMF). Its 
differential dF when the system evolves reversibly from u to u + du is the infinitesimal work that has to be provided 
to the system (or received from it) to allow such a reversible transformation. This work is that of the external forces 
added in the equations of motion of the local modeslS . 

The extension to the Gibbs free energy is straightforward^^: 

AG{N,{a},T;u) = G{N,{a},T;u)-Go{N,P,T) = - <fy2(Pr) _ -du' , (11) 

J — \ ' N.{a-},T-u' 
i 

in which the thermal averages are computed under fixed temperature, stress tensor and polarization. This yields: 



In this last expression, the relaxation of strain inducing a volume variation {fl{P, T, u)) along the path is taken into 
account. 



C. Computational details 

In this work, we have computed AG{N, {a} ,T;u) = G{N, {a} ,T;u) ~ Go{N, {a} , T) for BaTiOs, for 6 tempera- 
tures around the Curie temperature along the [100] direction. This is achieved through the following procedure. The 
calculation is performed every Am = 0.001 ao, starting from u = 0, which is found to be a fine enough grid to make 
the integration. Except in Sec. IIIII in which the effect of the supercell size is systematically tested, we use a 12 x 12 
X 12 supercell with periodic boundary conditions. 

First we perform for each point 10^ steps of MD by combining fixed pressure (Parrinello- Rahman algorithm), fixed 
temperature (Nose-Hoover algorithm) and fixed polarization^^. The (hydrostatic) pressure is fixed to -4.8 GPa, which 
is the value used by Zhong et al to correct the underestimation of the lattice constant of BTO within the Local Density 
Approximation^. The 50000 first steps are used to equilibrate the system. The last 50 000 are used to average the 
strain. The evolution of P computed at constant pressure and for a given temperature is indeed accompanied by an 
evolution of the strain, as can be seen on Fig. [2l 

Then, for each value of u of the grid, the calculation is restarted (for 10^ steps) under fixed strain tensor (the one 
averaged over the previous simulation), polarization and temperature. The final 50 000 steps of this second run are 
used to calculate the thermal average of the forces acting on the local modes. We check at the end of this second run 
that the stress tensor (obtained from the average over the 50 000 last steps) corresponds to the hydrostatic pressure 
of -4.8 GPa (it works with an accuracy of ~ 0.01 GPa). This procedure in two times is applied to ensure that the 
combination of the three algorithms (Nose-Hoover, Parrinello-Rahman, fixed polarization) does not induce artefacts 
that would bias the statistical averages. In fact, the free energy profiles obtained from the forces averaged over the 50 
000 last steps of the first run are almost identical to those obtained from the second one. We point out that at least 
50 000 steps are necessary to compute properly the mean force (forces are microscopic quantities having very large 
fluctuations). The mean force in numerically integrated as in Refs. [2oH25l . 



III. STATES OF CONSTRAINED POLARIZATION AT T=280 K (« T^ - 15-17 K) 

In the method employed here, the free energy is computed from MD simulations under fixed polarization, as 
explained above. The result is a sampling of phase space under the constraint of fixed order parameter, defined from 
an average over the whole supercell. However, at low order parameter, the system might separate into domains or at 
least become inhomogeneou&l£' inside the supercell. 

Moreover, the existence of domains would yield a difficulty to define properly the free energy as a volume quantity, 
by incorporating surface or interfacial contributions in the computed free energy. In Sec. lIII Al we check the reasonable 



7 



homogeneity of the states of constrained order parameter obtained in 12 x 12 x 12 supercells through 50000 time steps 
for the lowest temperature studied (280 K). In Sec. lIIIBl we test that this homogeneity is systematically improved by 
increasing the number of time steps of the simulation. The existence of interfacial contributions can be systematically 
tested by increasing the supercell size (Sec. lIIIC]) or decreasing the temperature (Sec. IIIID]) . 

A. Homogeneity of the macroscopic states of constrained order parameter 

We first perform a series of simulation tests in a 12 x 12 x 12 supercell: we check that all along the simulated 
path, the system does not separate into domains, even for small values of the order parameter. We compute the 
time-averaged local mode on each site of the simulation supercell: 

k=l 

in which Ui{k) is the local mode at cell i at time step k {N^t is the number of time steps), and plot their distribution 
(Fig. [1]). Clearly, a multidomain state would result in two distinct peaks at least for one component {x, y or z). In 
all cases, we obtain single-peak distributions. On the graphs of Fig. [H we have put arrows to indicate approximately 
the values of the spontaneous polarization of tetragonal BTO at the temperature considered here (T=280 K). A 
two-domain state would result in two peaks approximately localized close to those arrows. In all the cases, the width 
of the local mode distribution is well below this spontaneous polarization. 

The absence of ferroelectric domains is also obvious from the evolution of the strain tensor components along the 
simulated path, plotted on Fig. [2] for six temperatures around the Curie temperature: the evolution is smooth and 
continuous, with rji = r]2 = rj^ 0.012 for zero order parameter {u^ = Uy = Uz = 0). If ferroelectric domains 
were present, one component would jump to a high value corresponding roughly to the tetragonal strain at the 
corresponding temperature. Note that this 0.012 strain obtained at low order parameter is precisely the value that 
would be extrapolated from the paraelectric phase across Tc, suggesting that our constrained simulation for zero 
order parameter allows to reach some constrained paraelectric phase below the Curie temperature, which provides a 
physical picture fully compatible with Landau theory. 

Thus, the computed free energy does probably not contain significant interfacial contribution (we show hereafter in 
Sec, nil Cl that the present free energy does not significantly vary when increasing the supercell size to 16 x 16 x 16). 
Anyway, this has to be controlled severely in each simulation. A quick guess on what is going on can be obtained by 
examining the strain tensor components evolution (similarly, if one computes a Helmholtz energy under fixed strain, 
by examining the stress). 

Moreover, as clearly explained in Ref. [l^, if a multidomain state occurs at fixed low order parameter, the increase 
of the order parameter simply results in one domain growing at the expend of the other, or in a succession of more 
or less complex interface motions. As a consequence, as long as the domain walls move without interacting with each 
other, the free energy remains "fiat" and yields very characteristic curves^. Such phase separation might occur at 
low temperature and low displacivenessi^. Precisely, this does not occur in the present simulations (Fig. [TOl) . that are 
performed around the Curie temperature. 

Anyway such cases occur at "low" temperature, i.e. a few 10 K below the phase transition (probably when the 
correlation length falls below the supercell size - see hereafter), or when the geometry enforces ferroelectricity (for 
instance we find it in thin films or with strongly distorted cells, for which a strong electrostrictive coupling favours a 
ferroelectric state). 

B. Simulation time 

The simulations are systematically improved by increasing the simulation time: Fig.[3]shows that the time-averaged 
local mode distribution becomes more and more peaked as the simulation time is increased, showing that we are not 
sampling some metastable state. To emphasize this point, we plot on Fig.|4]the mean force per 5-atom cell (/i™) 
and the free energy curve (obtained by integrating the mean force) computed at T=280 K after averaging over 50 
000 steps and 500 000 steps. A systematic improvement is found. Any phase separation would result in plateaus 
in the mean force curve, which are not observed. Moreover, the average over 500 000 steps provides a very smooth 
and regular curve. We observe also that the potentials of mean force obtained in the two cases are almost identical, 
because the numerical errors on the forces tend to cancel when the integration is performed. Thus very good quality 
curves are obtained with an average over 50 000 steps. 
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FIG. 1: Distribution of the time-averaged components of the local modes (arb. units) (averaged over 50 000 steps), for 5 
constrained values of the order parameter, at T=280 K (panels a,b,c,d and e). The purple (resp. blue and orange) lines refer to 
the X (resp. y and z) component. The graphs show that the macroscopic states of constrained polarization remain homogeneous 
along the simulation (no phase separation). The broadness of the peaks reflect fluctuations and can be systematically reduced 
(see Fig. [S]). It can be compared to panel f, where the distribution of the local modes in an unconstrained system at T=310K 
{i.e. above Tc) is plotted after the same simulation time. The black arrows localize approximately the values of the spontaneous 
polarization at T=280 K (where the peaks are expected if phase separation occurs). The supercell is 12 x 12 x 12 and the 
local modes are in ao (lattice constants) units. 



C. EfFet of the supercell size 

The effect of the supercell size is tested by performing simulations at T=280 K in a 16 x 16 x 16 and in a 20 x 
20 X 20 supercell. The total mean force obtained as a function of mean local mode, as well as the corresponding 
potential of mean force (free energy) are shown on Fig. [H 

In the case of the 16 x 16 x 16 supercell, the agreement is very good with respect to the 12 x 12 x 12 cell: the 
mean force and free energy curves are very close to those presented above. Plotting the local mode distribution as 
well as the strain tensor component evolution indicates the absence of phase separation at low order parameter and 
confirms that interfacial contributions to the free energy are negligible, if any. 

However, in the 20 x 20 x 20 supercell, phase separation starts to manifest at low order parameter, yielding a 
characteristic trend on the mean force and free energy curve: some "plateaus" do appear in the mean force curve. This 
can be observed on the local mode distribution (Fig. [6]) as well as on the strain tensor components. Consequently, the 
mean force and the free energy deviate slightly from the two curves obtained with smaller supercells. More precisely, 
three regions can be distinguished in that case, separated on Fig. [5] by black arrows: at low order parameter (close to 
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FIG. 2: Diagonal components of homogeneous strain tensor 771 (red diamonds), 172 (orange diamonds) and rjs (blue diamonds) 
as a function of the order parameter, for six temperatures around Tc. P = -4.8 GPa, the Parrinello-Rahman algorithm is used 
to fix the stress along the path. 



zero), the mean force and the free energy are flat, indicating pseudo-domains with roughly opposite polarizations in 
the directions perpendicular to x {x is the direction of the fixed order parameter). Then at Ux ~ 0.005 ao, the mean 
force increases but through a linear curve (instead of a convex curve as in the smaller supercells) . From Fig. [6l this 
situation corresponds to pseudo-domains with polarization oriented along x, very inhomogeneous (Fig. El panels (f) 
and (g)), but the polarization in the perpendicular directions is maintained to its fixed value. Then at Ux ~ 0.015 ao, 
the distributions of the various components of the time-averaged local modes recover single-peak features, and the 
mean force recovers similar values as in smaller supercells. 

The free energy, as a consequence, even though the agreement with previous curves is still good (Fig. [71), tends to 
differ slightly. This trend is probably enstrengthed in larger supercells. 

Too large supercells thus lead to low order parameter states in which inhomogeneous configurations are dominating 
the sampHng, and the free energies obtained are consequently not smooth in these regions. In the following, we use 
the 12 X 12 X 12 supercell to compute the free energy around Tc. 

From the previous tests, we may assume that at T=280 K, the supercell sizes L=12 and L=16 are lower than 
some characteristic length of the system that can be compared to or interpreted as a correlation length. L=20 lattice 
constants is probably close to this correlation length. We show in the following that the computed free energies for 
L=12 are fully comparable to the phenomenological Landau free energies of the litterature. Supercell sizes ranging 
from L=12 to L=20 provide free energies (as a function of order parameter) that are quasi-similar at this temperature. 
Note that when approaching Tc from below, the correlation length increases and the free energy curves computed 
around Tc with L=12 are thus probably more sensitive to finite-size effects. 
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FIG. 3: Distribution of the time-averaged components of the local modes (arb. units), for 3 constrained values of the order 
parameter, at T=280 K. Panels a,b,c: average over 50000 steps. Panels d,e,f: average over 1000000 steps. The purple (resp. 
blue and orange) lines refer to the x (resp. y and z) component. The supercell is 12 x 12 x 12 the local modes are in ao 
(lattice constants) units. 



D. States of constrained polarization at lower temperature 

When the temperature decreases below 280 K, the states of constrained polarization at low order parameter become 
progressively more and more inhomogeneous, especially for T < 240 K. This is illustrated on Fig. [HI that represents 
the distribution of the time-averaged local modes in a L=12 supercell between T=270 K and T=210 K. 

At very low temperature and low order parameter, this distribution exhibits several peaks that reflect the appearance 
of very inhomogeneous states (Fig. [9]) . 



IV. RESULTS 



A. Computation of the free energy around Tc using a 12 X 12 X 12 supercell 

The Gibbs free energy curves as a function of u^^ computed in a 12 x 12 x 12 supercell for T=280, 290, 292.5, 
295, 297.5 and 300 K are shown on Fig. [TOl They correspond very well to the proflles expected from Landau theory 
for a first-order phase transition, with a free energy barrier separating the paraelectric phase from the ferroelectric 
phase for a few K below and above Tc. From these curves, we can localize the Curie temperature between 295 and 
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FIG. 4: Mean force as a function of mean local mode (arb. units) and potential of mean force (obtained by integration of minus 
the mean force), (a): averaging over 50 000 steps; (b): averaging over 500 000 steps. Note that the two free energy curves have 
not exactly the same depth due to linear deviations as the integration is performed. 



297.5 K (Iniguez et al localize it at 297 K with the same hamiltonian^) . The free energy barrier separating the cubic 
paraelectric phase and the tetragonal ferroelectric phase at Tc is AG ~ 0.012-0.015 meV/5-atom cell. The order of 
magnitude of this small barrier is in good agreement with that given by phenomenological Landau Gibbs free energies 
such as that of Wang et aU^ for which AG ~ 0.01 meV/5-atom cell at Tc. 



B. Comparison with phenomenological potentials 



We now compare the results with classical phenomenological Landau potentials, that assume in particular a Hnear 
dependence of the quadratic coefEcient with temperature. In the litterature, these potentials can be found as Gibbs 
free energies for the stress-free crystal, under the form of polynomial functions up to 6*'' or 8*'' order in the Px, Py 
and Pz- For example, the 8*'' order development writes: 
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FIG. 5: Mean force as a function of mean local mode (arb. units) and potential of mean force (obtained by integration of minus 
the mean force), (a): 16 x 16 x 16 supercell, average over 50 000 steps (b): 20 x 20 x 20 supercell, average over 50 000 steps. 



G(T; P) - a, {P^ + P^^ + P^) + an {P.^ + P,^ + P*) 
+ai2{P^P^ + P^Pl + P^Pj) + am {P', + + Pf) 
+an2{Pl{P^ + Pt) + Py{P' + Pt) + P^iP' + Py)) 
+ai23P^PSP! + aiiii(P| + P^ + P!) 
+ani2{P!{P^ + Pi) + PliPl + Pi) + P'APl + P^)) 

+an22{PtPt+PtPt + P'yPt) 

+an2^{P^P^Pl + P^P^Pl + PtP^P^) 

The free energy along the x axis {Px ^ 0, Py = Pz = 0) writes: 

G(T; P) = aiPl -f a^P^ + amP^ + annP^ (14) 

A fit of our free energy curves by G'^'-order or 8"*-order polynomial functions (without odd-order term) aP^ -f hP^ -f 
cP^{+dP^) should directly provide coefficients comparable to the ai, an, am and aim of the litterature. Thus 
we fit these curves by such polynomial functions and extract the coefficients (the volume variation along the curve is 
accounted for when performing the calculation). 
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FIG. 6: Distribution of the time-averaged components of the local modes (arb. units), for 4 constrained values of the order 
parameter, at T=280 K. Panels a,b,c,d: 16 x 16 x 16 supercell. Panels e,f,g,h: 20 x 20 x 20 supercell. The purple (resp. blue 
and orange) lines refer to the x (resp. y and z) component. The local modes are in ao (lattice constants) units. 



First the computed free energy curves fit very well by such functions. The quadratic term obtained by fitting is 
plotted on Fig. [11] as a function of T for the two fits (6*'' and S*'' order). On the same figure, we plot the quadratic 
term of the phenomenological potential of Wang et ape. ai = 3.61xlO^(T-To) (V.m.C-i). In Landau theory for first 
order phase transitions, Tq is, just below Tc, the temperature above which the paraelectric state {P — 0) becomes a 
(local) minimum of G. Wang et al set this value at 391 K, just below the real Curie temperature. In our case. To has 
to be readjusted below the Curie temperature of the Effective Hamiltonian. The best agreement with the S*'' order 
polynomial fit is obtained for To=275 K (green circles on Fig. [iT]), which is indeed just below the Curie temperature 
of the hamiltonian (« 297 K). For this value of Tq, the agreement between the phenomenological potential of Wang 
et a9^ and our work is very good. The temperature dependence of the quadratic term fitted is linear on the small 
range of temperature considered here (280-300 K). This would probably not be the case if the study was extended to 
lower temperatures^. 

The temperature evolutions of the other coefficients are shown on Figs.[l2] (quartic and sixth order coefficient) and 
[T3l (eight-order coefficient) in the case of a fit by an 8*'' order polynomial function. We find an excellent agreement 
close to Tc with the two phenomenological potentials of Wang et and Li et a^, for the three coefficients an, 
am and aim. 
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FIG. 7: Free energies as a function of mean local mode (meV/5-atom cell), as calculated for the three supercell sizes. Yellow 
diamonds: 12 x 12 x 12 (50 000 steps); Orange diamonds: 12 x 12 x 12 (500 000 steps); Red diamonds: 20 x 20 x 20 (50 
000 steps); Blue diamonds: 16 x 16 x 16 (50 000 steps). 



C. Natural decomposition of the free energy 

Since the effective hamiltonian decomposes into various contribution^ (in particular onsite, short-range, dipole- 
dipole, elastic- local mode interaction), the force acting on the z*'* local mode f}"^ also naturally decomposes into: 



dui 



jself _|_ j'short _|_ j'dpl _^ jint 



with obvious notations. This leads to a natural decomposition of the free energy Ai^(iV, {r/} , T; m) = 
F(7V,{77},T;tx)-Fo(^,M,T) as: 

Ai?(iV, {r;} ,T-u)^ AF^'^ (iV, {77} , T; u) + Ai?^''-* (TV, {,7} ,T-u) + AF'^p' [N, {77} , T; u) + Ai?"* {N, {77} , T; u), (15) 
which extends easily to the Gibbs free energy: 

AG{N,{a} ,T;P) = AG'^^^ (N, {a} ,T; P) + AG''""'\N, {a} ,T; P) + AG'^p' (N, {a} ,T- P) + AG"'\N, {a} ,T; P), 

(16) 
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FIG. 8: Distribution of the time-averaged local modes in a 12 x 12 x 12 supercell after 100 000 steps for a constrained value 
of the order parameter {ux — 0, Uy = 0, Uz = 0). The unit is the same as in the previous figures. The local modes are in ao 
(lattice constants) units. 



First, we notice that in the vicinity of the phase transition, these contributions are much higher (at least two 
orders of magnitude) than the free energy itself, confirming the well-known trend, according to which the existence of 
ferroelectricity is a very delicate balance between large destabilizing contributions (here short-range and onsite) and 
large stabiHzing ones (dipole-dipole and elastic- local mode interaction), as illustrated on Fig. [141 



1. Short-range and electrostatic parts 

Examining the various contributions as a function of u (mean local mode) , we notice something a priori surprising: 
among the four contributions, two of them (the short-range and the dipole-dipole) are independent on the temperature. 
This would not be the case if these contributions were expressed as a function of P instead of u (because the volume 
depends on u, see Fig. [21). In fact, this is simply due to the fact that the corresponding parts of the Effective 
Hamiltonian are quadratic in the local modes: 



ipshort \ ^ ja.f3 , rshort ^ \ ^ to./?,, 

^,3,a,0 J;I3 

Thus the forces vary linearly with the local modes, and their thermal average (in any ensemble) is 
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FIG. 9: Distribution of the time-averaged local modes in a 12 x 12 x 12 supercell after 350 000 steps at T=200 K for a 
constrained value of the order parameter {ux = 0, Uy = 0, itz =0). The unit is the same as in the previous figures. The local 
modes are in ao (lattice constants) units. 



Under fixed u, and under the hypothesis that the system is homogeneous (^j,f3, {uj^p) — up): 



Thus 
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E/ fshort\ ja.p 
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Up 



13 I i,3 



Let us introduce the 3x3 matrix: J^^p = Yl,i j '^ij^ ■ 



dG- 



short 



dUa 
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The same arguments stand of course in the canonical ensemble: 

Qp^short ^ 

It follows: 
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The SR contribution to the free energy is thus independent on the temperature and quadratic in the mean local 
mode. 
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For the dipole-dipole part, we have also: 

dQdpi^ Q pd-pi 

dUa dUa 

It follows: 



2X!Q«,/3"/3 (20) 

/3 



with Qa,f3 — J2i j Qtj^ Q matrix defined in Ref. l7|to compute the electrostatic energy. 

The previous equations allow to compute directly the short-range and dipole-dipole part of the free energy from 
the coefficients of the effective hamiltonian without simulation. These parts are quadratic in u and independent on 
the temperature. 

However, we stress that these two parts are independent on the temperature when expressed as a function of the 
mean local mode u. If they are expressed as a function of the polarization P, a temperature dependence appears, 
related to the variation of the cell volume {Vl/N) with temperature: P = N Z* /Q,{P,T,u)u. 



2. Onsite part and elastic-local mode interaction part 

The temperature dependence of the free energy, when expressed as a function of u, thus originates only in the 
onsite and elastic-local mode interaction parts of the Effective Hamiltonian. Let us examine the onsite energy: the 
corresponding part of the EflFective Hamiltonian is local and has an harmonic part and an anharmonic part: 

Prom the arguments given above, only the anharmonic part is likely to generate a temperature dependence of the 
corresponding part of the free energy: 



We have also: 
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With the notations of Ref. 0, we have: 

^F^^'f = Nn,\\uf + ^KnLrm (22) 

AG-'^=iV.^2Nf + AGrX™ (23) 

The contributions from anharmonic terms might contain also a priori an harmonic contribution to the free energy. 
Although this term in the Effective Hamiltonian only has quartic contributions, we have shown above that a fitting 
of the free energy on polynomial functions may contain higher-order terms. 

Finally, we focus on the elastic-local mode interaction energy: 

Bi^a,l3mi^)Ut,aU^^i3, 

with r]i{i) — rj^ + ?yf(i) (homogeneous strain + inhomogeneous strain). The same treatment can not be applied 
here due to presence of the inhomogeneous strain. This contribution is thus a priori dependent on the temperature 
and this is observed in our simulations. 

V. DISCUSSION AND CONCLUSION 

In this work, we have computed the free energy of barium titanate as a function of polarization directly from 
an effective hamiltonian. We have used Molecular Dynamics combined with Thermodynamic Integration. This is 
equivalent to define the free energy as the potential of the mean force acting on the polarization. One obtains, for 
a given temperature, the difference of free energy between P = and P. The simulations can be performed under 
fixed volume/strain or fixed pressure/stress conditions, which gives access to Helmholtz free energies or to Gibbs 
free energies. The gradient of the free energy is related to the thermal average of the total forces, computed under 
the constraint of fixed polarization. From the decomposition of the forces into the different contributions coming 
from the Effective Hamiltonian, a natural decomposition of the free energy of a ferroelectric into four parts (onsite, 
dipole-dipole, short-range and elastic-local mode interaction) has been suggested. In particular we have shown that 
the dipole-dipole and short-range contributions to the free energy are, in the framework of the Effective Hamiltonian 
of Zhong et af , independent on the temperature when expressed as a function of u instead of P and are quadratic in 
u. 

The computation has been performed for a set of temperatures around Tc and by using a 12 x 12 x 12 supercell. 
Although the comparison of such a free energy with a Landau potential is still controversial™, a very good agreement 
has been found with Landau phenomenological potentials of the litteratur o^^i^^ aroud Tc. We have checked that in 
the range of temperature considered, phase separation at low order parameter does not occur within this supercell. 
Moreover we have shown, by systematically increasing the supercell size and the number of time steps of the simulation, 
that the computed free energy does not contain significant interfacial contributions and can be considered as a volume 
(extensive) quantity. We think that these techniques can be useful to approach the Landau free energy of a ferroelectric 
near the Curie temperature and by using atomic-scale simulations. 

As reminded in the introduction, the concept of Landau free energy requires the choice of a spatial averaging length 
L to define the order parameter as a smooth and continuous physical quantity. If this averaging length is chosen 
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above the correlation length ^, some local modes might be unsufficiently correlated within the averaging region and 
microscopic configurations with phase separation are likely to appear frequently in the sampling when a low value of 
the (averaged) order parameter is fixed. These phenomena, already reported in the framework of the (j)^ model, have 
been also observed in the present case at low temperature or when the supercell size is very large. Such considerations 
lead to the concept of coarse-grained free energiesi^ and elegant approaches have been proposed to circumvent the 
corresponding difficulties^. 

However, as reminded also, the definition of the correlation volume is rather complex in ferroelectric systems due 
to the very anisotropic nature of the correlations. Further investigation of the notion of correlation length/ volume in 
ferroelectrics should be of primary interest to enlight the concept of Landau free energy in those materials. 
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the thermodynamic integration providing the Gibbs free energy is formally identical to that providing the Helmholtz free 

energy. This is due to the peculiar form of the Effective Hamiltonian, in which the local modes and the homogeneous strain 

are independent variables. 



